Adj=[];
corrcoef=1.3;
%% adjustible step size with Non-negative concentrations control
% Reminder: nRsteps = 100;      % take nRsteps reaction steps per diffusion step
% REMINDER: dt2 = dt/nRsteps;   % reaction time-step
% tStart3  = tic;
react_t=0;
adj_count=0;
while (react_t<=dt)
    dts_i=dt2;
    [Cdc42Th, BemGEF42h, BemGEFh, BemGEFch, Cdc42Dh, GDI42h, GDI42ch, GDIch, Cdc42Dich, GDI42ich, Recich,Rech,RecAh] = ...
         Euler_step_GDI42ic_mixedCyt(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,Recic,Rec,RecA,    dts_i);

   % Species_h  = Species    + dts_i*react_flux(dts,Species);
    while (~isempty(find(Cdc42Th<0))  | ~isempty(find(BemGEF42h<0)) | ~isempty(find(BemGEFh<0)) ...
          | ~isempty(find(BemGEFch<0)) | ~isempty(find(Cdc42Dh<0))   | ...
            ~isempty(find(GDI42h<0))   | ~isempty(find(GDI42ch<0))   | ~isempty(find(GDIch<0)) ...
          | ~isempty(find(GDI42ich<0)) | ~isempty(find(Cdc42Dich<0)) | ~isempty(find(Rech<0)) ...
          | ~isempty(find(RecAh<0)) | ~isempty(find(Recich<0)) )
        dts_i=dts_i/1.5;
        adj_count=adj_count+1; %Masha debug
        [Cdc42Th, BemGEF42h, BemGEFh, BemGEFch, Cdc42Dh, GDI42h, GDI42ch, GDIch, Cdc42Dich, GDI42ich, Recich,Rech,RecAh] = ...
            Euler_step_GDI42ic_mixedCyt(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,Recic,Rec,RecA,    dts_i);
    end     
    react_t=react_t+dts_i;
    Cdc42T   = Cdc42Th;
    BemGEF42 = BemGEF42h;
    BemGEF   = BemGEFh;
    BemGEFc  = BemGEFch;
    Cdc42D   = Cdc42Dh;
    GDI42    = GDI42h;
    GDI42c   = GDI42ch;
    GDIc     = GDIch;
    Cdc42Dic = Cdc42Dich;
    GDI42ic  = GDI42ich;
    Recic= Recich;
    Rec = Rech;
    RecA     = RecAh;
end 
Adj=[Adj adj_count];
%tElapsed3=[tElapsed3 toc(tStart3)];

%for nnt = 1:nRsteps
%%NO step adjustment, simple one step integration
    %     [Cdc42T, BemGEF42, BemGEF, BemGEFc, Cdc42D, GDI42, GDI42c, GDIc, Cdc42Dic, GDI42ic, Recic,Rec,RecA] = ...
    %              Euler_step_GDI42ic_mixedCyt(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,Recic,Rec,RecA,    dt2);    
%% Adjustable step size with error estimation
%      lp = 0;
%      react_t = 0;
%      et = 0.0001;
%      adj_count=0; %Masha debug
%      
%      while  lp ~= 2
%          adj_count=adj_count+1; %Masha debug
%          %Masha dev sept 2012 Added Recic,Recn,RecAn species to all Euler function calls 
%         [Cdc42Th, BemGEF42h, BemGEFh, BemGEFch, Cdc42Dh, GDI42h, GDI42ch, GDIch, Cdc42Dich, GDI42ich, Recich,Rech,RecAh] = ...
%             Euler_step_GDI42ic_mixedCyt(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,Recic,Rec,RecA,    dt2);
%         [Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2, Cdc42Dich2,GDI42ich2,Recich2,Rech2,RecAh2] = ...
%             Euler_step_GDI42ic_mixedCyt(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,Recic,Rec,RecA,    dt2/2);
%         [Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2, Cdc42Dich2,GDI42ich2,Recich2,Rech2,RecAh2] = ...
%             Euler_step_GDI42ic_mixedCyt(Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2,Cdc42Dich2,GDI42ich2,Recich2,Rech2,RecAh2,dt2/2);
%         
%         Cdc42Te    = max(max(abs(Cdc42Th   - Cdc42Th2   )));
%         BemGEF42e  = max(max(abs(BemGEF42h - BemGEF42h2 )));
%         BemGEFe    = max(max(abs(BemGEFh   - BemGEFh2   )));
%         BemGEFce   = max(max(abs(BemGEFch  - BemGEFch2  )));
%         Cdc42De    = max(max(abs(Cdc42Dh   - Cdc42Dh2   )));
%         GDI42e     = max(max(abs(GDI42h    - GDI42h2    )));
%         GDI42ce    = max(max(abs(GDI42ch   - GDI42ch2   )));
%         GDIce      = max(max(abs(GDIch     - GDIch2     )));
%         Cdc42Dice  =        abs(Cdc42Dich  - Cdc42Dich2 )  ;
%         GDI42ice   =        abs(GDI42ich   - GDI42ich2  )  ;
%         Recice=        abs(Recich-Recich2)  ; %Masha dev sept2012
%         Rece  =max(max(abs(Rech  -  Rech2))); %Masha dev sept2012
%         RecAe      =max(max(abs(RecAh      -      RecAh2))); %Masha dev sept2012
%         
%         r = max([Cdc42Te,BemGEF42e,BemGEFe,BemGEFce,Cdc42De,GDI42e,GDI42ce,GDIce,Cdc42Dice,GDI42ice,Recice,Rece,RecAe]);
% 
%         % doing the step again with different stepsize        
%         if r > et && lp ~= 2
%             dt2 = dt2/corrcoef;     
% 
%             if lp == 1
%                 lp = 0;
%             end
%             % moving on
%         else
% 
%             Cdc42T   = Cdc42Th2;
%             BemGEF42 = BemGEF42h2;
%             BemGEF   = BemGEFh2;
%             BemGEFc  = BemGEFch2;
%             Cdc42D   = Cdc42Dh2;
%             GDI42    = GDI42h2;
%             GDI42c   = GDI42ch2;
%             GDIc     = GDIch2;
%             Cdc42Dic = Cdc42Dich2;
%             GDI42ic  = GDI42ich2;
%             Recic= Recich2;
%             Rec = Rech2;
%             RecA     = RecAh2;
%             
%             % adding the step on 
%             react_t = react_t + dt2;
%             % setting new stepsize for next step
%             dt2 = corrcoef*dt2;
%             
%             if lp == 1
%                 lp = 2;
%             end
%             
%         end
%             
%         % Making sure we stay within the dt interval
%         if react_t + dt2 > dt && lp ~= 2
%             dt2 = dt - react_t;
%             lp = 1;
%         end            
%             
%         
%         %checking for negative values and reducing tolerance 
%         if lp == 2
%             if (~isempty(find(Cdc42T<0))  | ~isempty(find(BemGEF42<0)) | ~isempty(find(BemGEF<0)) ...
%              | ~isempty(find(BemGEFc<0)) | ~isempty(find(Cdc42D<0))   | ...
%                ~isempty(find(GDI42<0))   | ~isempty(find(GDI42c<0))   | ~isempty(find(GDIc<0)) ...
%              | ~isempty(find(GDI42ic<0)) | ~isempty(find(Cdc42Dic<0)) | ~isempty(find(Rec<0)) ...
%              | ~isempty(find(RecA<0)) )
%                 lp = 0;
%                 react_t = 0;
%                 et = et/10;
%                 fprintf(neg_store,' %5.4f \t%5.4e \t%5.4e \n', curr_t,et*10,et);  
%             end
%         end
%      end
%      %display([' ' num2str(adj_count)]); %Masha debug
%      Adj = [Adj adj_count];
%end
 %display(['Number of dt2 adjustments ' num2str(mean(Adj)) ' +- ' num2str(std(Adj))]);
    %% diffusion in cargo step